The present invention relates to power-flow/voltage control in utility/industrial power networks of the types including many power plants/generators interconnected through transmission/distribution lines to other loads and motors. Each of these components of the power network is protected against unhealthy or alternatively faulty, over/under voltage, and/or over loaded damaging operating conditions. Such a protection is automatic and operates without the consent of power network operator, and takes an unhealthy component out of service by disconnecting it from the network. The time domain of operation of the protection is of the order of milliseconds.
The purpose of a utility/industrial power network is to meet the electricity demands of its various consumers 24-hours a day, 7-days a week while maintaining the quality of electricity supply. The quality of electricity supply means the consumer demands be met at specified voltage and frequency levels without over loaded, under/over voltage operation of any of the power network components. The operation of a power network is different at different times due to changing consumer demands and development of any faulty/contingency situation. In other words healthy operating power network is constantly subjected to small and large disturbances. These disturbances could be consumer/operator initiated, or initiated by overload and under/over voltage alleviating functions collectively referred to as security control functions and various optimization functions such as economic operation and minimization of losses, or caused by a fault/contingency incident.
For example, a power network is operating healthy and meeting quality electricity needs of its consumers. A fault occurs on a line or a transformer or a generator which faulty component gets isolated from the rest of the healthy network by virtue of the automatic operation of its protection. Such a disturbance would cause a change in the pattern of power flows in the network, which can cause over loading of one or more of the other components and/or over/under voltage at one or more nodes in the rest of the network. This in turn can isolate one or more other components out of service by virtue of the operation of associated protection, which disturbance can trigger chain reaction disintegrating the power network.
Therefore, the most basic and integral part of all other functions including optimizations in power network operation and control is security control. Security control means controlling power flows so that no component of the network is over loaded and controlling voltages such that there is no over voltage or under voltage at any of the nodes in the network following a disturbance small or large. As is well known, controlling electric power flows include both controlling real power flows which is given in MWs, and controlling reactive power flows which is given in MVARs. Security control functions or alternatively overloads alleviation and over/under voltage alleviation functions can be realized through one or combination of more controls in the network. These involve control of power flow over tie line connecting other utility network, turbine steam/water/gas input control to control real power generated by each generator, load shedding function curtails load demands of consumers, excitation controls reactive power generated by individual generator which essentially controls generator terminal voltage, transformer taps control connected node voltage, switching in/out in capacitor/reactor banks controls reactive power at the connected node.
Control of an electrical power system involving power-flow control and voltage control commonly is performed according to a process shown in FIG. 4. The various steps entail the following.    Step-10: Obtain on-line/simulated readings of open/close status of all switches and circuit breakers, and read data of maximum and minimum reactive power generation capability limits of PV-node generators, maximum and minimum tap positions limits of tap changing transformers, and maximum power carrying capability limits of transmission lines, transformers in the power network, or alternatively read data of operating limits of power network components.    Step-20: Obtain on-line readings of real and reactive power assignments/schedules/specifications/settings at PQ-nodes, real power and voltage magnitude assignments/schedules/specifications/settings at PV-nodes and transformer turns ratios. These assigned/set values are controllable and are called controlled variables/parameters.    Step-30: Resulting voltage magnitudes and angles at power network nodes, power flows through various power network components, reactive power generations by generators and turns ratios of transformers in the power network are determined by performance of loadflow computation, for the assigned/set/given/known values from step-20 or previous process cycle step-60, of controlled variables/parameters.    Step-40: The results of Loadflow computation of step-30 are evaluated for any over loaded power network components like transmission lines and transformers, and over/under voltages at different nodes in the power system    Step-50: If the system state is acceptable implying no over loaded transmission lines and transformers and no over/under voltages, the process branches to step-70, and if otherwise, then to step-60    Step-60: Changes the controlled variables/parameters set in step-20 or as later set by the previous process cycle step-60 and returns to step-30    Step-70: Actually implements the corrected controlled variables/parameters to obtain secure/correct/acceptable operation of power system
Overload and under/over voltage alleviation functions produce changes in controlled variables/parameters in step-60 of FIG. 5. In other words controlled variables/parameters are assigned or changed to the new values in step-60. This correction in controlled variables/parameters could be even optimized in case of simulation of all possible imaginable disturbances including outage of a line and loss of generation for corrective action stored and made readily available for acting upon in case the simulated disturbance actually occurs in the power network. In fact simulation of all possible imaginable disturbances is the modern practice because corrective actions need be taken before the operation of individual protection of the power network components.
It is obvious that loadflow computation consequently is performed many times in real-time operation and control environment and, therefore, efficient and high-speed loadflow computation is necessary to provide corrective control in the changing power system conditions including an outage or failure of any of the power network components. Moreover, the loadflow computation must be highly reliable to yield converged solution under a wide range of system operating conditions and network parameters. Failure to yield converged loadflow solution creates blind spot as to what exactly could be happening in the network leading to potentially damaging operational and control decisions/actions in capital-intensive power utilities.
The power system control process shown in FIG. 5 is very general and elaborate. It includes control of power-flows through network components and voltage control at network nodes. However, the control of voltage magnitude at connected nodes within reactive power generation capabilities of electrical machines including generators, synchronous motors, and capacitor/inductor banks, and within operating ranges of transformer taps is normally integral part of loadflow computation as described in “LTC Transformers and MVAR violations in the Fast Decoupled Load Flow, IEEE Trans., PAS-101, No. 9, PP. 3328-3332, September 1982.” If under/over voltage still exists in the results of loadflow computation, other control actions, manual or automatic, may be taken in step-60 in the above and in FIG. 5. For example, under voltage can be alleviated by shedding some of the load connected.
The prior art and present invention are described using the following symbols and terms:
Ypq=Gpq'jBpq: (p-q) th element of nodal admittance matrix without shunts
Ypp=Gpp+jBpp: p-th diagonal element of nodal admittance matrix without shunts
yp=gp+jbp: total shunt admittance at any node-p
Vp=ep+jfp=Vp<θp: complex voltage of any node-p
Δθp, ΔVp: voltage angle, magnitude corrections
Δep, Δfp: real, imaginary components of voltage corrections
Pp+jQP: net nodal injected power calculated
ΔPP+jΔQp: nodal power residue or mismatch
RPP+jRQp: modified nodal power residue or mismatch
Φp: rotation or transformation angle
[RP]: vector of modified real power residue or mismatch
[RQ]: vector of modified reactive power residue or mismatch
[Yθ]: gain matrix of the P-θ loadflow sub-problem defined by eqn. (1)
[YV]: gain matrix of the Q-V loadflow sub-problem defined by eqn. (2)
m: number of PQ-nodes
k: number of PV-nodes
n=m+k+1: total number of nodes
q>p: q is the node adjacent to node-p excluding the case of q=p
[ ]: indicates enclosed variable symbol to be a vector or a matrix
PQ-node: load-node, where, Real-Power-P and Reactive-Power-Q are specified
PV-node: generator-node, where, Real-Power-P and Voltage-Magnitude-V are specified Bold lettered symbols represent complex quantities in description.
    Loadflow Computation: Each node in a power network is associated with four electrical quantities, which are voltage magnitude, voltage angle, real power, and reactive power. The loadflow computation involves calculation/determination of two unknown electrical quantities for other two given/specified/scheduled/set/known electrical quantities for each node. In other words the loadflow computation involves determination of unknown quantities in dependence on the given/specified/scheduled/set/known electrical quantities.    Loadflow Model: a set of equations describing the physical power network and its operation for the purpose of loadflow computation. The term ‘loadflow model’ can be alternatively referred to as ‘model of the power network for loadflow computation’. The process of writing Mathematical equations that describe physical power network and its operation is called Mathematical Modeling. If the equations do not describe/represent the power network and its operation accurately the model is inaccurate, and the iterative loadflow computation method could be slow and unreliable in yielding converged loadflow computation. There could be variety of Loadflow Models depending on organization of set of equations describing the physical power network and its operation, including Decoupled Loadflow Models, Super Decoupled Loadflow (SDL) Models, Fast Super Decoupled Loadflow (FSDL) Model, and Novel Fast Super Decoupled Loadflow (NFSDL) Model.    Loadflow Method: sequence of steps used to solve a set of equations describing the physical power network and its operation for the purpose of loadflow computation is called Loadflow Method, which term can alternatively be referred to as ‘loadflow computation method’ or ‘method of loadflow computation’. One word for a set of equations describing the physical power network and its operation is: Model. In other words, sequence of steps used to solve a Loadflow Model is a Loadflow Method. The loadflow method involves definition/formation of a loadflow model and its solution. There could be variety of Loadflow Methods depending on a loadflow model and iterative scheme used to solve the model including Decoupled Loadflow Methods, Super Decoupled Loadflow (SDL) Methods, Fast Super Decoupled Loadflow (FSDL) Method, and Novel Fast Super Decoupled Loadflow (NFSDL) Method. All decoupled loadflow methods described in this application use either successive (1θ, 1V) iteration scheme or simultaneous (1V, 1θ), defined in the following.
Prior art methods of loadflow computation of the kind carried out as step-30 in FIG. 5, are well known Gauss-Seidel Loadflow (GSL) and Super Super Decoupled Loadflow (SSDL) methods. The GSL method is well known to be not able to converge to high accuracy solution because of its iteration scheme that lacks self iterations, which realization led to the invention of Gauss-Seidel-Patel Loadflow (GSPL) method. The prior art methods will now the described.
Gauss-Seidel Loadflow: GSL
The complex power injected into the node-p of a power network is given by the following relation,
                                                        P              p                        -                          j              ⁢                                                          ⁢                              Q                p                                              =                                                    V                p                *                            ⁢                                                ∑                                      q                    =                    1                                    n                                ⁢                                                      Y                    pq                                    ⁢                                      V                    q                                                                        =                                                            V                  p                  *                                ⁢                                  Y                  pp                                ⁢                                  V                  p                                            +                                                V                  p                  *                                ⁢                                                      ∑                                          q                      >                      p                                                        ⁢                                                            Y                      pq                                        ⁢                                          V                      q                                                                                                          ⁢                                  ⁢                  Where          ,                                    (        1        )                                          P          p                =                  Re          ⁢                                          ⁢                      {                                          V                p                *                            ⁢                                                ∑                                      q                    =                    1                                    n                                ⁢                                                      Y                    pq                                    ⁢                                      V                    q                                                                        }                                              (        2        )                                          Q          p                =                              -            Im                    ⁢                                          ⁢                      {                                          V                p                *                            ⁢                                                ∑                                      q                    =                    1                                    n                                ⁢                                                      Y                    pq                                    ⁢                                      V                    q                                                                        }                                              (        3        )            Where, words “Re” means “real part of” and words “Im” means “imaginary part of”.
The Gauss-Seidel (GS) method is used to solve a set of simultaneous algebraic equations iteratively. The GSL-method calculates complex node voltage from any node-p relation (1) as given in relation (4).
                              V          p                =                              [                                          {                                                      (                                                                  PSH                        p                                            -                                              j                        ⁢                                                                                                  ⁢                                                  QSH                          p                                                                                      )                                    /                                      V                    p                    *                                                  }                            -                                                ∑                                      q                    >                    p                                                  ⁢                                                      Y                    pq                                    ⁢                                      V                    q                                                                        ]                    /                      Y            pp                                              (        4        )            Iteration Process
Iterations start with the experienced/reasonable/logical guess for the solution. The reference node also referred to as the slack-node voltage being specified, starting voltage guess is made for the remaining (n−1)-nodes in n-node network. Node voltage value is immediately updated with its newly calculated value in the iteration process in which one node voltage is calculated at a time using latest updated other node voltage values. A node voltage value calculation at a time process is iterated over (n−1)-nodes in an n-node network, the reference/slack node voltage being specified not required to be calculated. Now, for the iteration-(r+1), the complex voltage calculation at node-p equation (4) and reactive power calculation at node-p equation (3), becomes
                              V          p                      (                          r              +              1                        )                          =                              [                                                                                                      {                                                                        (                                                                                    PSH                              p                                                        -                                                          j                              ⁢                                                                                                                          ⁢                                                              QSH                                p                                                                                                              )                                                /                                                                              (                                                          V                              p                              *                                                        )                                                    r                                                                    }                                        -                                                                                                                                                                  ∑                                                  q                          =                          1                                                                          p                          -                          1                                                                    ⁢                                                                        Y                          pq                                                ⁢                                                  V                          q                                                      (                                                          r                              +                              1                                                        )                                                                                                                -                                                                  ∑                                                  q                          =                                                      p                            +                            1                                                                          n                                            ⁢                                                                        Y                          pq                                                ⁢                                                  V                          q                          r                                                                                                                                          ]                    /                      Y            pp                                              (        5        )                                          Q          p                      (                          r              +              1                        )                          =                              -            Im                    ⁢                      {                                                                                (                                          V                      p                      *                                        )                                    r                                ⁢                                                      ∑                                          q                      =                      1                                                              p                      -                      1                                                        ⁢                                                            Y                      pq                                        ⁢                                          V                      q                                              (                                                  r                          +                          1                                                )                                                                                                        -                                                                    (                                          V                      q                      *                                        )                                    r                                ⁢                                                      ∑                                          q                      =                      p                                        n                                    ⁢                                                            Y                      pq                                        ⁢                                          V                      q                      r                                                                                            }                                              (        6        )            Convergence
The iteration process is carried out until changes in the real and imaginary parts of the set of (n−1)-node voltages calculated in two consecutive iterations are all less than the specified tolerance—ε, as shown in relations (7) and (8). The lower the value of the specified tolerance for convergence check, the greater the solution accuracy.|ΔeP(r+1)|=|eP(r+1)−ePr|<ε  (7)|ΔfP(r+1)|=|fP(r+1)−fPr|<ε  (8)Accelerated Convergence
The GS-method being inherently slow to converge, it is characterized by the use of an acceleration factor applied to the difference in calculated node voltage between two consecutive iterations to speed-up the iterative solution process. The accelerated value of node-p voltage at iteration-(r+1) is given byVp(r+1)(accelerated)=Vpr+β(Vp(r+1)−VPr)  (9)Where, β is the real number called acceleration factor, the value of which for the best possible convergence for any given network can be determined by trial solutions. The GS-method is very sensitive to the choice of β, causing very slow convergence and even divergence for the wrong choice.Scheduled or Specified Voltage at a PV-Node
Of the four variables, real power PSHp and voltage magnitude VSHp are specified at a PV-node. If the reactive power calculated using VSHp at the PV-node is within the upper and lower generation capability limits of a PV-node generator, it is capable of holding the specified voltage at its terminal. Therefore the complex voltage calculated by relation (5) by using actually calculated reactive power Qp in place of QSHp is adjusted to specified voltage magnitude by relation (10).Vp(r+1)=(VSHpVp(r+1))/|Vp(r+1)|  (10)Calculation Steps of Gauss-Seidel Loadflow (GSL) Method
The steps of loadflow computation method, GSL method are shown in the flowchart of FIG. 1a. Referring to the flowchart of FIG. 1a, different steps are elaborated in steps marked with similar numbers in the following. The words “Read system data” in Step-1 correspond to step-10 and step-20 in FIG. 5, and step-14, step-20, step-32, step-44, step-50 in FIG. 6. All other steps in the following correspond to step-30 in FIG. 5, and step-60, step-62, and step-64 in FIG. 6.    1. Read system data and assign an initial approximate solution. If better solution estimate is not available, set specified voltage magnitude at PV-nodes, 1.0 p.u. voltage magnitude at PQ-nodes, and all the node angles equal to that of the slack-node angle, which is referred to as the flat-start.    2. Form nodal admittance matrix, and Initialize iteration count r=1    3. Scan all the node of a network, except the slack-node whose voltage having been specified need not be calculated. Initialize node count p=1, and initialize maximum change in real and imaginary parts of node voltage variables DEMX=0.0 and DFMX=0.0    4. Test for the type of a node at a time. For the slack-node go to step-12, for a PQ-node go to the step-9, and for a PV-node follow the next step.    5. Compute Qp(r+1) for use as an imaginary part in determining complex schedule power at a PV-node from relation (6) after adjusting its complex voltage for specified value by relation (10)    6. If Qp(r+1) is greater than the upper reactive power generation capability limit of the PV-node generator, set QSHp=the upper limit Qpmax for use in relation (5), and go to step-9. If not, follow the next step.    7. If Qp(r+1) is less than the lower reactive power generation capability limit of the PV-node generator, set QSHp=the lower limit Qpmin for use in relation (5), and go to step-9. If not, follow the next step.    8. Compute Vp(r+1) from relation (5) using QSHp=Qp(r+1), and adjust its value for specified voltage at the PV-node by relation (10), and go to step-10    9. Compute Vp(r+1) from relation (5)    10. Compute changes in the real and imaginary parts of the node-p voltage by using relations (7) and (8), and replace current value of DEMX and DFMX respectively in case any of them is larger.    11. Calculate accelerated value of Vp(r+1) by using relation (9), and update voltage by Vpr=Vp(r+1) for immediate use in the next node voltage calculation.    12. Check for if the total numbers of nodes-n are scanned. That is if p<n, increment p=p+1, and go to step-4. Otherwise follow the next step.    13. If DEMX and DFMX both are not less than the convergence tolerance (ε) specified for the purpose of the accuracy of the solution, advance iteration count r=r+1 and go to step-3, otherwise follow the next step    14. From calculated and known values of complex voltage at different power network nodes, and tap position of tap changing transformers, calculate power flows through power network components, and reactive power generation at PV-nodes.Decoupled Loadflow
In a class of decoupled loadflow methods, each decoupled method comprises a system of equations (11) and (12) differing in the definition of elements of [RP], [RQ], and [Yθ] and [YV]. It is a system of equations for the separate calculation of voltage angle and voltage magnitude corrections.[RP]=[Yθ][Δθ]  (11)[RQ]=[YV][ΔV]  (12)Successive (1θ, 1V) Iteration Scheme
In this scheme (11) and (12) are solved alternately with intermediate updating. Each iteration involves one calculation of [RP] and [Δθ] to update [θ] and then one calculation of [RQ] and [ΔV] to update [V]. The sequence of relations (13) to (16) depicts the scheme.[Δθ]=[Yθ]−1[RP]  (13)[θ]=[θ]+[Δθ]  (14)[ΔV]=[YV]−1[RQ]  (15)[V]=[V]+[ΔV]  (16)
The scheme involves solution of system of equations (11) and (12) in an iterative manner depicted in the sequence of relations (13) to (16). This scheme requires mismatch calculation for each half-iteration; because [RP] and [RQ] are calculated always using the most recent voltage values and it is block Gauss-Seidel approach. The scheme is block successive, which imparts increased stability to the solution process. This in turn improves convergence and increases the reliability of obtaining solution.
Super Super Decoupled Loadflow: SSDL
This method is not very sensitive to the restriction applied to nodal transformation angles; SSDL restricts transformation angles to the maximum of −48 degrees determined experimentally for the best possible convergence from non linearity considerations, which is depicted by relations (19) and (20). However, it gives closely similar performance over wide range of restriction applied to the transformation angles say, from −36 to −90 degrees.RPp=(ΔPp Cos Φp+ΔQp Sin Φp)/Vp2 —for PQ-nodes  (17)RQp=(ΔQp Cos Φp−ΔPp Sin Φp)/Vp —for PQ-nodes  (18)Cos Φp=Absolute(Bpp/v(Gpp2+Bpp2))≧Cos(−48°)  (19)Sin Φp=−Absolute(Gpp/v(Gpp2+Bpp2))≧Sin(−48°)  (20)RPp=ΔPp/(KpVp2) —for PV-nodes  (21)Kp=Absolute(Bpp/Yθpp)  (22)
                              Y          ⁢                                          ⁢                      θ            pq                          =                  [                                                                                          -                                    ⁢                                      Y                    pq                                                                                                                                          -                                        ⁢                    for                    ⁢                                                                                  ⁢                    branch                    ⁢                                                                                                              ⁢                                                                                                            ⁢                                          r                      /                      x                                        ⁢                                                                                  ⁢                    ratio                                    ≤                  3.0                                                                                                                          -                                    ⁢                                      (                                                                  B                        pq                                            +                                              0.9                        ⁢                                                  (                                                                                    Y                              pq                                                        -                                                          B                              pq                                                                                )                                                                                      )                                                                                                                                          -                                        ⁢                    for                    ⁢                                                                                  ⁢                    branch                    ⁢                                                                                  ⁢                                          r                      /                      x                                        ⁢                                                                                                              ⁢                                                                                                            ⁢                    ratio                                    >                  3.0                                                                                                                          -                                    ⁢                                      B                    pq                                                                                                                                                                                                                        -                                                    ⁢                          for                          ⁢                                                                                                          ⁢                          branches                          ⁢                                                                                                          ⁢                          connected                          ⁢                                                                                                          ⁢                          between                          ⁢                                                                                                          ⁢                          two                                                                                                                                                              PV                          ⁢                                                      -                                                    ⁢                          nodes                          ⁢                                                                                                          ⁢                          or                          ⁢                                                                                                                                            ⁢                                                                                                                                          ⁢                          a                          ⁢                                                                                                          ⁢                          PV                          ⁢                                                      -                                                    ⁢                          node                          ⁢                                                                                                          ⁢                          and                          ⁢                                                                                                          ⁢                          the                          ⁢                                                                                                          ⁢                          slack                          ⁢                                                      -                                                    ⁢                          node                                                                                                      ⁢                                                                                                                                                (        23        )                                          YV          pq                =                  [                                                                                          -                                    ⁢                                      Y                    pq                                                                                                                                          -                                        ⁢                    for                    ⁢                                                                                  ⁢                    branch                    ⁢                                                                                                              ⁢                                                                                                            ⁢                                          r                      /                      x                                        ⁢                                                                                  ⁢                    ratio                                    ≤                  3.0                                                                                                                          -                                    ⁢                                      (                                                                  B                        pq                                            +                                              0.9                        ⁢                                                  (                                                                                    Y                              pq                                                        -                                                          B                              pq                                                                                )                                                                                      )                                                                                                                                          -                                        ⁢                    for                    ⁢                                                                                  ⁢                    branch                    ⁢                                                                                  ⁢                                          r                      /                      x                                        ⁢                                                                                                              ⁢                                                                                                            ⁢                    ratio                                    >                  3.0                                                                                        (        24        )                                          Y          ⁢                                          ⁢                      θ            pp                          =                                            ∑                              q                >                p                                      ⁢                                          -                Y                            ⁢                                                          ⁢                              θ                pq                            ⁢                                                          ⁢              and              ⁢                                                          ⁢                              YV                pp                                              =                                    b              p              ′                        +                                          ∑                                  q                  >                  p                                            ⁢                              -                                  YV                  pq                                                                                        (        25        )                                                                                    b                p                ′                            =                                                (                                                                                    QSH                        p                                            ⁢                      Cos                      ⁢                                                                                          ⁢                                              Φ                        p                                                              -                                                                  PSH                        p                                            ⁢                      Sin                      ⁢                                                                                          ⁢                                                                        Φ                          p                                                /                                                  V                          s                          2                                                                                                      )                                -                                                      b                    p                                    ⁢                  Cos                  ⁢                                                                          ⁢                                      Φ                    p                                    ⁢                                                                          ⁢                  or                                                                                                                        b                p                ′                            =                              2                ⁢                                                      (                                                                                            QSH                          p                                                ⁢                        Cos                        ⁢                                                                                                  ⁢                                                  Φ                          p                                                                    -                                                                        PSH                          p                                                ⁢                        Sin                        ⁢                                                                                                  ⁢                                                  Φ                          p                                                                                      )                                    /                                      V                    s                    2                                                                                                          (        26        )            bp′=(QSHp Cos Φp−PSHp Sin Φp/Vs2)−bp Cos Φp orbp′=2(QSHp Cos Φp−PSHp Sin Φp)/Vs2  (26)where, Kp is defined in relation (22), which is initially restricted to the minimum value of 0.75 determined experimentally; however its restriction is lowered to the minimum value of 0.6 when its average over all less than 1.0 values at PV nodes is less than 0.6. Restrictions to the factor Kp as stated in the above is system independent. However it can be tuned for the best possible convergence for any given system. In case of systems of only PQ-nodes and without any PV-nodes, equations (23) and (24) simply be taken as Yθpq=YVpq=−Ypq.
Branch admittance magnitude in (23) and (24) is of the same algebraic sign as its susceptance. Elements of the two gain matrices differ in that diagonal elements of [YV] additionally contain the b′ values given by relations (25) and (26) and in respect of elements corresponding to branches connected between two PV-nodes or a PV-node and the slack-node. Relations (19) and (20) with inequality sign implies that transformation angles are restricted to maximum of −48 degrees for SSDL. The method consists of relations (11) to (26). In two simple variations of the SSDL method, one is to make YVpq=Yθpq and the other is to make Yθpq=YVpq.
Calculation Steps of Super Super Decoupled Loadflow (SSDL) Method
The steps of loadflow computation method, SSDL method are shown in the flowchart of FIG. 1b. Referring to the flowchart of FIG. 1b, different steps are elaborated in steps marked with similar letters in the following. The words “Read system data” in Step-1 correspond to step-10 and step-20 in FIG. 5, and step-14, step-20, step-32, step-44, step-50 in FIG. 6. All other steps in the following correspond to step-30 in FIG. 5, and step-60, step-62, and step-64 in FIG. 6.    a. Read system data and assign an initial approximate solution. If better solution estimate is not available, set voltage magnitude and angle of all nodes equal to those of the slack-node. This is referred to as the slack-start.    b. Form nodal admittance matrix, and Initialize iteration count ITRP=ITRQ=r=0    c. Compute Cosine and Sine of nodal rotation angles using relations (19) and (20), and store them. If they, respectively, are less than the Cosine and Sine of −48 degrees, equate them, respectively, to those of −48 degrees.    d. Form (m+k)×(m+k) size matrices [Yθ] and [YV] of (11) and (12) respectively each in a compact storage exploiting sparsity. The matrices are formed using relations (23), (24), (25), and (26). In [YV] matrix, replace diagonal elements corresponding to PV-nodes by very large value (say, 10.0**10). In case [YV] is of dimension (m×m), this is not required to be performed. Factorize [Yθ] and [YV] using the same ordering of nodes regardless of node-types and store them using the same indexing and addressing information. In case [YV] is of dimension (m×m), it is factorized using different ordering than that of [Yθ].    e. Compute residues [ΔP] at PQ- and PV-nodes and [ΔQ] at only PQ-nodes. If all are less than the tolerance (ε), proceed to step-n. Otherwise follow the next step.    f. Compute the vector of modified residues [RP] as in (17) for PQ-nodes, and using (21) and (22) for PV-nodes.    g. Solve (13) for [Δθ] and update voltage angles using, [θ]=[θ]+[Δθ].    h. Set voltage magnitudes of PV-nodes equal to the specified values, and Increment the iteration count ITRP=ITRP+1 and r=(ITRP+ITRQ)/2.    i. Compute residues [ΔP] at PQ- and PV-nodes and [ΔQ] at PQ-nodes only. If all are less than the tolerance (ε), proceed to step-n. Otherwise follow the next step.    j. Compute the vector of modified residues [RQ] as in (18) for only PQ-nodes.    k. Solve (15) for [ΔV] and update PQ-node magnitudes using [V]=[V]+[ΔV]. While solving equation (15), skip all the rows and columns corresponding to PV-nodes.    l. Calculate reactive power generation at PV-nodes and tap positions of tap changing transformers. If the maximum and minimum reactive power generation capability and transformer tap position limits are violated, implement the violated physical limits and adjust the loadflow solution.    m. Increment the iteration count ITRQ=ITRQ+1 and r=(ITRP+ITRQ)/2, and Proceed to step-e.    n. From calculated and known values of voltage magnitude and voltage angle at different power network nodes, and tap position of tap changing transformers, calculate power flows through power network components, and reactive power generation at PV-nodes.